6 - Working with Rasters¶

In this notebook, we will:

  • Download rasters from EarthByte's webDAV server
  • Plot rasters
  • Resize and respace rasters
  • Linearly interpolate point data on rasters.

Import all needed packages, and create PlateReconstruction and Plot objects for the Muller et al. (2019) plate tectonic model.

In [1]:
import os

import cartopy.crs as ccrs
import gplately
import matplotlib.pyplot as plt
import numpy as np

Using GPlately's download module, we can easily download a variety of plate tectonic models!

In [2]:
gdownload = gplately.download.DataServer("Muller2019")
rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()
coastlines, continents, COBs = gdownload.get_topology_geometries()

time = 0  # Ma
model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
gplot = gplately.PlotTopologies(model, time, coastlines, continents, COBs)

Checking whether the requested files need to be updated...
Requested files are up-to-date!
Checking whether the requested files need to be updated...
Requested files are up-to-date!

GPlately's DataServer for rasters¶

Let's use GPlately's DataServer object to download Muller et al. 2019 netCDF age grids.

There is a unique age grid for each millionth year - let's access the 0 Ma age grid by passing time to get_age_grid. It is returned as a masked array.

In [3]:
time = 0  # Ma
muller_2019_age_grid = gdownload.get_age_grid(time)

Checking whether the requested files need to be updated...
Requested files are up-to-date!

Plotting rasters¶

The Raster object allows us to work with age grids and other rasters. Let's specify that muller_2019_age_grid is an array and visualise the data with imshow.

In [4]:
graster = gplately.Raster(model, data=muller_2019_age_grid, extent=[-180, 180, -90, 90])

fig = plt.figure(figsize=(16, 12))
plt.imshow(graster.data, origin="lower", cmap="YlGnBu")
plt.show()

Let's plot this netCDF grid along with coastlines, mid-ocean ridges and subduction zones (with teeth) resolved from the Muller et al. 2019 plate model.

In [5]:
fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
gplot.time = time
gplot.plot_grid(ax1, graster.data, cmap='YlGnBu', vmin=0, vmax=200)
gplot.plot_coastlines(ax1, edgecolor='k', facecolor='1', alpha=0.1)
gplot.plot_trenches(ax1, color='magenta', zorder=5)
gplot.plot_ridges(ax1, color='b', zorder=5)
gplot.plot_transforms(ax1, color='lightblue', linewidth=0.75)
gplot.plot_subduction_teeth(ax1, color='magenta')
ax1.set_title("{} Ma".format(time))
plt.show()

Resizing and resampling rasters¶

Let's resize and resample the the present-day Muller et al. 2019 agegrid. We can do this using resize

In [6]:
model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
muller_2019_age_grid = gdownload.get_age_grid(time)
graster = gplately.Raster(model, data=muller_2019_age_grid, extent=[-180, 180, -90, 90])

# Set grid size in x and y directions
graster.resize(20, 10, overwrite=True)

# Plot resampled and resized age grid
fig = plt.figure(figsize=(16,12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude = 20))
gplot.plot_grid(ax1, graster.data, cmap='YlGnBu', vmin=0, vmax=200)
gplot.plot_coastlines(ax1, edgecolor='k', facecolor='1', alpha=0.1)
gplot.plot_trenches(ax1, color='magenta', zorder=5)
gplot.plot_ridges(ax1, color='b', zorder=5)
gplot.plot_subduction_teeth(ax1, color='magenta')
ax1.set_title("{} Ma".format(time))
plt.show()

Checking whether the requested files need to be updated...
Requested files are up-to-date!

Downloading general rasters¶

Let's visualise an ETOPO1 relief raster. GPlately also makes it easy to import, using get_raster

Alternatively, you can import a local netcdf file using gplately.grids.read_netcdf_grid

In [7]:
etopo = gdownload.get_raster("ETOPO1_tif")
print(np.shape(etopo))
print(np.size(etopo) // 3)

Checking whether the requested files need to be updated...
Requested files are up-to-date!
(2700, 5400, 3)
14580000

etopo is a large (5400 x 2700 pixels) RGB image, with a total of 14,580,000 grid points. We can visualise it using a number of methods, including Raster.imshow.

In [8]:
raster = gplately.Raster(data=etopo, origin="upper")

fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
raster.imshow(ax=ax1, interpolation="none")
plt.show()

We can also resize this RGB raster using the Raster.resize and Raster.resample methods. Here we resample it to a $2\degree \times 2\degree$ resolution, then plot it using Raster.imshow.

In [9]:
etopo_downscaled = gplately.Raster(
    data=raster.resample(2, 2),
    origin="upper",
)
etopo_downscaled.imshow(projection=ccrs.Robinson(central_longitude=180), interpolation="none")
Out[9]:
<matplotlib.image.AxesImage at 0x2d76b1430>

Reconstructing rasters¶

The ETOPO1 raster can be reconstructed back in time using the Raster and PlateReconstruction classes.

First we create our Raster object, where we can specify our plate model we want to reconstruct with.

After this, we use reconstruct to reconstruct the raster to a given time.

In [10]:
graster = gplately.Raster(
    plate_reconstruction=model,
    data=etopo,
    extent="global",  # equivalent to (-180, 180, -90, 90)
    origin="upper",  # or set extent to (-180, 180, 90, -90)
)
white_rgb = (255, 255, 255)  # RGB code for white, to fill gaps in output
rgb = graster.reconstruct(50, threads=4, fill_value=white_rgb)

fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=180))
gplot.plot_grid(ax1, rgb, origin="upper")
Out[10]:
<matplotlib.image.AxesImage at 0x34827a460>

Rasters can be reconstructed in-place (inplace=True), and fill_value can be set to any valid matplotlib colour when reconstructing RGB images.

In [11]:
etopo_flipped = np.flipud(etopo)

raster = gplately.Raster(
    plate_reconstruction=model,
    data=etopo_flipped,
    extent="global",
    origin="lower",
)

raster.reconstruct(time=75, threads=4, fill_value="darkblue", inplace=True)
raster.imshow(projection=ccrs.Robinson())
Out[11]:
<matplotlib.image.AxesImage at 0x34362ea90>

By default, [Raster.reconstruct](reconstruct uses self.plate_reconstruction.static_polygons to assign plate IDs to grid points. To override this behaviour, pass any collection of pygplates.Feature (e.g. list, pygplates.FeatureCollection, etc) to the partitioning_features argument.

In [12]:
raster = gplately.Raster(
    plate_reconstruction=model,
    data=etopo,
    extent="global",
    origin="upper",
)

raster.reconstruct(140, partitioning_features=continents, threads=4, fill_value="grey", inplace=True)
raster.imshow(projection=ccrs.Orthographic(0, -80))
plt.gca().set_title("Reconstructed to 140 Ma")
Out[12]:
Text(0.5, 1.0, 'Reconstructed to 140 Ma')

Rasters can be also be reverse reconstructed forward in time!

In [13]:
raster = gplately.Raster(model, data=etopo, time=0, origin="upper")

reconstructed = raster.reconstruct(50, fill_value="white", threads=4)

raster_reconstructed = gplately.Raster(model, data=reconstructed, time=50, origin="upper")

reverse_reconstructed = raster_reconstructed.reconstruct(0, fill_value="white", threads=4)

# plot
fig, axs = plt.subplots(3, 1, figsize=(8, 12), subplot_kw={"projection": ccrs.Mollweide(central_longitude=0)})
axs[0].imshow(etopo, transform=ccrs.PlateCarree())
axs[1].imshow(reconstructed, transform=ccrs.PlateCarree())
axs[2].imshow(reverse_reconstructed, transform=ccrs.PlateCarree())

axs[0].set_title("Original image (present day)")
axs[1].set_title("Reconstructed to 50 Ma")
axs[2].set_title("Reverse reconstructed back to present day")
plt.show()

Reconstructing netCDF rasters¶

Similar to above, we can also reconstruct netcdf grids! GPlately also includes ETOPO1 as a netcdf for download.

In [14]:
etopo = gdownload.get_raster("ETOPO1_grd").astype("float")

graster = gplately.Raster(
    plate_reconstruction=model, 
    data=etopo,
    extent="global",
    origin="lower",  # note: here it is 'lower', while the ETOPO tiff uses 'upper'
)
reconstructed_etopo = graster.reconstruct(50, threads=4)

# plot
fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
gplot.plot_grid(ax1, reconstructed_etopo, origin="lower")
plt.show()

Checking whether the requested files need to be updated...
Requested files are up-to-date!

The reconstructed netcdf can also be saved using gplately.grids.write_netcdf_grid

In [15]:
# Save the reconstructed ETOPO grid to netCDF 
# (note: this ETOPO netCDF grid is high-res, ~190+MB)
save_filename = os.path.join(
    "NotebookFiles",
    "reconstructed_etopo.nc",
)
gplately.grids.write_netcdf_grid(save_filename, reconstructed_etopo)

Linear interpolation¶

As an example of linear interpolation, let's project points along trenches out a distance of 250 km and check which projected points lie inside a continental raster.

In [16]:
# Define the continental grid
continental_grid_filename = os.path.join(
    "NotebookFiles",
    "continental_grid_0.nc",
)
continental_raster = gplately.Raster(
    model,
    continental_grid_filename,
    extent=[-180, 180, -90, 90],
)

# Extract subduction polarity angle, and the lat-lon coordinates of each trench point 
trench_data = model.tesselate_subduction_zones(time, ignore_warnings=True)
trench_normal_azimuthal_angle = trench_data[:, 7]
trench_pt_lon = trench_data[:, 0]
trench_pt_lat = trench_data[:, 1]
    
# Project the 250km arc-trench distance onto an Earth sphere
arc_distance = 250 / (gplately.tools.geocentric_radius(trench_pt_lat) / 1e3)
    
# Lat and lon coordinates of all trench points after being projected out 250 km in the direction of subduction.
dlon = arc_distance * np.sin(np.radians(trench_normal_azimuthal_angle))
dlat = arc_distance * np.cos(np.radians(trench_normal_azimuthal_angle))
ilon = trench_pt_lon + np.degrees(dlon)
ilat = trench_pt_lat + np.degrees(dlat)

Let's use GPlately to linearly interpolate these projected trench points onto the continental grids using Raster.interpolate

In [17]:
# Use GPlately to interpolate these new points on the defined grid
sampled_points = continental_raster.interpolate(ilon, ilat, method='linear')

# sampled_point[0] is a list of points in the grid (ascribed the integer 1). Collect their indices. 
in_raster_indices = sampled_points > 0

# Get the lat-lon coordinates of the in_raster points
lat_in = ilat[in_raster_indices]
lon_in = ilon[in_raster_indices]

Plot the in-raster points along with the raster, trenches, and coastlines.

In [18]:
fig = plt.figure(figsize=(16, 12), dpi=100)
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
gplot.time = time
gplot.plot_grid_from_netCDF(ax1, continental_grid_filename, cmap="twilight", alpha=0.5, vmin=0, vmax=200)
gplot.plot_coastlines(ax1, edgecolor='k', facecolor='1', alpha=0.1)
gplot.plot_trenches(ax1, color='r', zorder=5)

# Plot the trench points in-continent
ax1.scatter(lon_in, lat_in, s=0.25, transform=ccrs.PlateCarree(), label="Points in continental grids")
plt.title("%i Ma" %time)
plt.legend(loc="lower left")
plt.show()